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ABSTRACT 

This paper reports the progress made towards 
developing complete blood flow simulation capability in human, 
especially, in the presence of artificial devices such as valves 
and ventricular assist devices. Devices modeling poses unique 
challenges different from computing the blood flow in natural 
hearts and arteries. There are many elements needed to quantify 
the flow in these devices such as flow solvers, geometry 
modeling including flexible walls, moving boundary procedures 
and physiological characterization of blood. As a first step, 
computational technology developed for aerospace applications 
was extended to the analysis and development of a ventricular 
assist device (VAD), i.e. a blood pump. The blood flow in a VAD 
is practically incompressible and Newtonian, and thus an 
incompressible Navier-Stokes solution procedure can be applied. 
A primitive variable formulation is used in conjunction with the 
overset grid approach to handle complex moving geometry. The 
primary purpose of developing the incompressible flow analysis 
capability was to quantify the flow in advanced turbopump for 
space propulsion system. The same procedure has been 
extended to the development of NASA-DeBakey VAD that is 
based on an axial blood pump. Due to massive computing 
requirements, high-end computing is necessary for simulating 
three-dimensional flow in these pumps. Computational, 
experimental and clinical results are presented. 

INTRODUCTION 

Approximately 20 million people worldwide suffer 
annually from congestive heart failure (CHF), a quarter of them in 
America alone. In the United States, an alarmingly low 2,000 to 
2,500 donor hearts are available each year. One potential 


approach to improve this situation is to use a mechanical device 
to boost or to create blood flow in patients suffering from 
hemodynamic deterioration; that is, loss of blood pressure and 
lowered cardiac output. The goal of this device can be to replace 
the natural heart, i.e. total artificial heart, or to assist an ailing 
heart, i.e. ventricular assist device (VAD) [1]. In either approach, 
the device can be used to bridge the gap while waiting for a 
matching donor heart for transplantation. However, to ease the 
shortage of donor hearts, making these devices suitable for long- 
term or permanent use would be an ultimate goal. 

Another benefit of an assist device is the potential for 
providing time for the natural heart to recover. In some patients, 
it has been observed that the natural heart can recover by 
unloading the pumping requirement through the use of a VAD. 

In what conditions this might happen is not very well quantified 
at this time and should involve physiological particulars of 
patients among other factors. The challenge is to design a device 
that can deliver the required blood circulation while not 
adversely impacting human physiological conditions. In this 
paper, since computational aspects for developing either total 
artificial heart or an assist device are largely the same, mechanical 
blood pumping devices are represented generically by a VAD. 

Requirements for a VAD related to fluid dynamics are 
demanding such as; simplicity and reliability; small size for ease 
of implantation; pumping capacity to supply 5 liter/min of blood 
against 100 mmHg pressure; high pumping efficiency to minimize 
power requirements; minimum hemolysis and thrombus 
formation. In addition to fluid dynamic issues, there are many 
other important aspects to be taken care of such as material 
compatibility to human, controls and implantation procedures. 
Due to complexity of flow physics and delicate operating 
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conditions, empirical approach to quantify the flow phenomena 
in a VAD is not all that straightforward, very time consuming and 
expensive, especially to study many design variations. Here 
comes in the potential utility of computational simulation tools 
for the development of these devices. In this paper, the 
discussion is focused on how fluid dynamic issues of VAD can 
be resolved via computational approach. 

Computational flow analysis of blood flow through 
mechanical devices such as VADs and artificial valves is 
extremely challenging [2]. Flow is unsteady and involves moving 
parts. For a complete analysis of a VAD such as the axial blood 
pump illustrated in Figure 1, human circulatory system simulation 
(e.g. Quarteroni [3]) has to be coupled to the device in use. 
However, for the purpose of developing mechanical components, 
a truncated circulation system can be modeled. For example, 
empirical inflow condition can be specified at the inlet of a VAD. 
Even with this type of simplifications, computational approach 
can produce flow field data in great detail, thus shedding lights 
to obtain a better understanding of the dominant flow physics 
produced by an artificial device. Especially, computational 
analysis can be utilized to optimize the design of mechanical 
devices at a significantly lower cost and time than required by an 
empirical approach. In addition to the geometric and operational 
complexities, these devices introduce a variety of flow 
phenomena, which are not normally existing in a natural heart. 
These include transition, turbulence, boundary layer separation, 
rotational effects, tip vortex and reverse flow phenomena. 

In the present study, a computational procedure for 
developing a VAD is discussed. Since the computational tools 
are derived from those developed under rocket propulsion 
systems development, tools development on rocket turbopump 
will be discussed first followed by a presentation of applications 
to VAD flows. 


FORMULATION AND METHOD OF SOLUTION 

As illustrated in Figure 1, an assist device is placed to 
draw blood from an ailing heart and then to pump the blood into 
the aorta. Blood flow in large vessels as in this arrangement is 
believed to behave as Newtonian fluid. The flow can thus be 
solved by the following incompressible Navier-Stokes equations: 


the relative reference frame is moving around the x-axis, the 
source term, S, is given by 
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Where Q is the rotational speed. Relative velocity components 
are written in terms of absolute velocity comp onents u a v a w a . 
u =u a 

v = v a + Qz 
w = w a - Qy 

For component analysis, rotational steady formulation may 
provide quick yet valuable approximate solutions. Solution 
procedures and a family of flow solver codes, INS3D, has been 
developed at NASA Ames Research Center [4-6]. This code has 
been applied to numerous viscous internal flow problems 
including turbopump component analysis for liquid rocket 
engine. 



Figure 1. Schematic of a VAD based on an axial-flow pump. 
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Where t is the time, the Cartesian coordinates, U[ the 

corresponding velocity components, P the pressure, and ij the 
viscous-stress tensor. When the equations are solved in a 
steady rotating frame of reference, the centrifugal and Coriolis 
forces are added as source terms to the governing equations. If 


The algorithm used in the current work is based on the 
method of pseudo-compressibility where the time -derivative of 
pressure is introduced into the continuity equation. The elliptic - 
parabolic type partial differential equations are transformed into 
the hyperbolic -parabolic type. Since the convective terms of the 
resulting equations are hyperbolic, upwind differencing can be 
applied to these terms. Third - and fifth-order flux-difference 
splitting is used for the convective terms. The upwind 
differencing leads to a more diagonally dominant system than 
does central differencing and does not require any user specified 
artificial dissipation. The viscous flux derivatives are computed 
using central differencing. In the steady-state formulation, the 
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time derivatives are differenced using the Euler backward 
formula. The equations are solved iteratively in pseudo-time until 
the solution converges to a steady state. In the time -accurate 
formulation, the equations are iterated to convergence in 
pseudo-time for each physical time step until the divergence of 
the velocity field is reduced below a specified tolerance value. 
Discretization of the governing equations results in a matrix 
equation that is solved iteratively using a non-factored Gauss- 
Seidel type line relaxation scheme that allows the use of a large 
pseudo -time step. Details of the numerical method can be found 
in the references cited. 

CFD TECHNOLOGY FOR LIQUID ROCKET PUMP 

Until recently, the high performance pump design 
process was not significantly different from that of 30 years ago. 
During the past 30 years a vast amount of experimental and 
operational experience has demonstrated that there are many 
important features of pump flows, which are not accounted for in 
the semi -empirical design process. Pumps being designed today 
are no mo're technologically advanced than those designed for 
the Space Shuttle Main Engine (SSME). During that same time 
span huge strides have been made in computers, in numerical 
algorithms, and in physical modeling. 

Rocket pumps involve full and partial blades, tip leakage and exit 
boundary to diffuser. In addition to the geometric complexities, a 
variety of flow phenomena are encountered in turbopump flows. 
These include turbulent boundary layer separation, wakes, 
transition, tip vortex resolution, three-dimensional effects, and 
Reynolds numbers effects. In order to increase the role of 
Computational Fluid Dynamics (CFD) in the design process, the 
CFD analysis tools must be evaluated and validated so that 
designers gain confidence in their use. 

As a part of extending the CFD technology to pump design 
work, INS3D code has been validated for pump component 
analysis using data from a rocket pump inducer [7]. The resulting 
computational procedure was applied to the flow through the 
SSME High Pressure Fuel Turbopump impeller and to the 
development of an advanced pump impeller (Kiris and Kwak, 

1994 [8]). The results from the advanced pump impeller flow 
analysis are presented next. 

In Figure 2, a cross sectional view of an advanced impeller is 
shown schematically. The computational model of an advanced 
pump includes the impeller and the exit shroud cavity region. 
Figure 3 shows the computational grid near the hub region of the 
impeller. The impeller design flow rate is 1,205 gal/min with a 
design speed of 6,322 rpm. The Reynolds Number for this 
calculation was 1 8 1,283 per inch. In Figure 4, the meridional 
velocity is shown at the impeller discharge. 


| 
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Figure 2. Schematic of an advanced pump impeller cross-section. 

A relative x-distance is measured from the shroud to hub, where 
x=1.0 is the hub. The meridional velocities, Cm, were integrated 
along a radial strip for each constant x- position and they were 
non-dimensionalized by the wheel speed of 249.5 ft/sec. 



Figure 3. Advanced pump impeller computational grid on the 
hub surface. 

The meridional velocity distribution for 5% and 10% recirculation 
from the exit shroud cavity were also plotted. When the exit 
shroud cavity has leakage to the impeller eye, the velocity peak 
at the impeller exit moves toward the center of the b2 width, 
where b2 is defined as the blade height at the impeller exit (see 
figure 2). However, the shroud leakage has only minor effects on 
the solution at r/rtip= 1.0275 (Figure 4). 
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In Figure 4, the symbols represent experimental data, and the 
lines represent Cm distributions for the flow with vaneless space 
at the exit of the impeller. 



Relative X at Impeller Exit 


Figure 4. Comparison of circumferentially averaged meridional 
velocity at the impeller exit. 

The test data shows that the peak is closer to the center of the b2 
width. The discrepancy between the computed results and 
experimental data is partially due to the recirculation flow in the 
hub cavity. The leakage at the hub cavity leads to a stronger 
recirculation region, which shifts the velocity peak to the center 
of the b2 width. Since the CFD analysis did not include the 
leakage at the hub cavity, the predicted recirculation region in 
the vaneless space is not as strong as in the experimental study. 

To simulate the entire turbopump including inducer, 
impeller, diffuser and volute, it is necessary to compute fully 
three-dimensional unsteady flow with components in relative 
motion. In the present study, the unsteady computing procedure 
for the entire turbopump has been developed using overset grid 
systems [9]. This capability is intended to provide a 
computational framework for analyzing rocket engine fuel and 
oxidizer supply system. Three-dimensional unsteady flow 
analysis with multiples of moving and stationary components 
require many physical time steps as well as a large number of grid 
points to resolve flow features of interest. Therefore, accelerating 
the computational procedure is of significant importance in 
engineering practice. One possibility of speed up comes from 
enhancements in computer hardware platforms, which at this time 
relies on parallel architecture. Advances in algorithms and 
efficient parallel implementations contribute to the other part of 
the speed-up. In the following, key elements in unsteady pump 
flow simulation are presented. 

The pump used in liquid propulsion system may include 
various rotating and stationary components where the flow can 
be axial or centrifugal and extremely unsteady. A stand-alone 
inducer is often used to validate the flow solver used. In Figure 
6, computed result of an inducer alone at a constant rotational 


speed is shown. A similar inducer was used later in VAD 
development. 



Surface Pressure 


Figure 6. Surface pressure for a pump inducer. 



Figure 7. Geometry of a SSME pump including inlet guide 
vanes, impeller and diffuser blades. 



Figure 8. An overset grid system for the impeller blade section 
with tip clearance. 

When rotating and stationary parts are included as shown in 
Figure 7, unsteady interactions between stationary and rotating 
components need to be included. To resolve the 
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complex geometry in relative motion, an overset grid approach is 
employed in the current effort. In this approach, a geometrically 
complex body is decomposed into a number of simple grid 
components. To illustrate this feature, the overset grid set up for 
the region around an impeller blade is shown in Figure 8. In order 
to solve the entire configuration including inlet guide vanes, 
impeller blades and diffuser blades, an overset grid system was 
constructed using 34.3 Million grid points with 1 14 zones in the 
present study. Since neighboring grids may overlap arbitrarily, 
these grids can be created independently from each other. Thus 
it is easier to generate high quality and nearly orthogonal grids 
using this approach. Connectivity between neighboring grids is 
established by interpolation at the grid outer boundaries. 
Addition of new components to the system and simulation of 
arbitrary relative motion between multiple bodies are achieved by 
establishing new connectivity without disturbing the existing 
grids. 

Parallel computing strategies vary depending on 
computer architecture such as memory arrangement relative to 
processing units. Two approaches have been implemented in the 
present study: the first approach is hybrid MPI/OpenMP and the 
second one is Multi Level Parallelism (MLP) developed at 
NASA-Ames Research Center. The first approach is obtained by 
using massage-passing interface (MPI) for inter-zone parallelism, 
and by using OpenMP directives for intra-zone parallelism. 
INS3D-MPIis based on the explicit massage-passing interface 
across MPI groups and is designed for coarse grain parallelism. 
The primary strategy is to distribute the zones across a set of 
processors. During the iteration, all the processors would 
exchange boundary data between processors whose zones 
shared interfaces with zones on other processors. A simple 
master-worker architecture was selected because it is relatively 
simple to implement and it is a common architecture for parallel 
CFD applications. All I/O was performed by master MPI process 
and data was distributed to the workers. After the initialization 
phase is complete, the program begins its main iteration loop. 

The second approach is obtained by using NAS-MLP 
routines. This approach differs from the MPI/OpenMP approach 
in a fundamental way in that it does not use messaging at all. All 
data communication at the coarsest and finest level is 
accomplished via direct memory referencing instructions. This 
approach also provides a simpler mechanism than MPI for 
converting legacy code, such as INS3D. For shared memory 
MLP, the coarsest level parallelism is supplied by spawning of 
independent processes via the standard UNIX fork. The 
advantage of the UNIX fork over MPI procedure is that the user 
does not have to change the initialization section of the large 
production code. Library of routines is used to initiate forks, to 
establish shared memory arenas, and to provide synchronization 
primitives. The boundary data for the overset grid system is 
archived in the shared memory arena by each process. Other 
processes access the data from the arena as needed. Figure 9 
shows the scalability for SSME impeller computations using 19.2 
Million grid points. 



Figure 9. Time (sec) per iteration for SSME impeller 
computations. 

Using both MPI/OpenMP hybrid approach and NAS-MLP 
parallel implementation, time -accurate computations of a pump 
shown in Figure 7 have been carried out on SGI Origin 2000 and 
3000 platforms. Initially the flow was at rest and the impeller 
started to rotate impulsively. An instantaneous snapshot of 
particle traces and pressure surfaces from these computations 
are shown in Figure 10 after three full impeller rotations were 
completed. This procedure provides flow-field details not readily 
available from experiments. A full-length version of this work will 
be presented in a separate paper. 



Figure 10. A snapshot of particle traces and pressure surfaces 
from unsteady turbopump computations. 

BLOOD PUMP AS A VENTRICULAR ASIST DEVICE 

In 1989, DeBakey Heart Center of the Baylor College of 
Medicine (BCM) began developing a new implantable VAD 
system jointly with NASA Johnson Space Center (JSC) [1]. This 
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LVAD is based on a fast rotating axial pump. In order to deliver 
the required blood flow rate, the rotational speed of this VAD is 
in the range of rocket pump operating condition. Therefore, the 
CFD procedures described earlier in conjunction with the rocket 
pump simulation was applied to this VAD development. By 
probing through computed results, regions of critical design 
interest can be identified such as regions of high turbulence 
shear flow which can damage the red blood cells and regions of 
re -circulation where blood clots may be formed. The ability to 
predict these phenomena expedited the development of the VAD 

First, the baseline VAD impeller was analyzed by 
solving the incompressible Navier-Stokes equations in a steady 
rotating frame of reference. Zonal multiblock grids were used in 
this component analysis. 



Figure 11. Computational grid for the baseline model of DeBakey 
VAD. 

As shown in Figure 1 1, the computational domain is divided into 
five zones with grid dimensions of 127 x 39 x 33, 127 x 39 x 33, 59 x 
21 x 7, 47 x 21 x 5, and 59 x 21 x 7, respectively. Zone 1 is the 
region between the suction side of the partial blade and the 
pressure side of the full blade; the region between the pressure 
side of the partial blade and the suction side of the full blade is 
filled by zone 2; and zones 3 through 5 allow tip -leakage effects 
to be included in the computational study and occupy the 
regions between the impeller blade tip and the casing. At the 
zonal interfaces, grid points were matched one-to-one. For all 
zones, an H-H type grid topology was used. An H-type surface 
grid was generated for each surface using an elliptic grid 
generator. The interior region of the three-dimensional grid was 
filled using an algebraic grid generator coupled with an elliptic 
smoother. Periodic boundary conditions were used at the end 
points in the rotational direction. The design flow of this impeller 
is 5 liters per minute and the design speed is 12,600 revolutions 
per minute (rpm). The problem was non-dimensionalized by the 
tube diameter (0.472 inches) and the impeller tip velocity. The 
solution was considered converged when the maximum residual 
had dropped at least five orders of magnitude. 

A parametric study was performed to optimize the 
impeller blade shape and the tip clearance. The desgn shown in 
Figure 1 1 was analyzed with different tip clearances. However, 
the performance of the baseline design could not be brought up 
to the level required for human implantation both in blood cell 
damage level and the thrombus formation. 


A modified designed was developed adding an inducer 
similar to the ones used in rocket pump. In collaboration with 
Micromed Technologies, NASA -JSC and BCM researchers, a 
new design consisting of the baseline impeller plus an inducer 
was investigated. Detail flow features and pump performances 
are compared for the two designs. The pressure gradient across 
the blades, and the pressure rise from inflow to outflow are 
compared in Figure 12. Blade geometry was then optimized in 
conjunction with detailed flow analysis. 



Figure 12. Pressure surfaces of the baseline design (top) and 
new impeller design (bottom). 



Figure 13. Meridional velocity distribution along impeller blade 
height of various designs. 

In Figure 13, the circumferentially averaged meridional velocity 
distribution is shown along the blade height for various designs. 
The original blade design is referred to as Design I. Design II has 
less blade curvature than Design I in the trailing edge region, and 
Design III has more blade curvature than Design I. In Design IV, 
the blade shape for Design I is kept and the tip clearance is 
reduced. In Design V, the hub region has the blade shape for 
Design I and the tip region has the blade shape for Design II. In 
this design, the impeller blades have backward lean near the 
trailing edge region. In Design VI, the blades have forward lean 
which includes Design III in the hub region and Design I in the 
tip region. Design VII has small tip clearance gap, Design I blade 
shape and an inducer geometry upstream of impeller blades. In 
Figure 13, all designs except Design VII showed back flow near 
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the hub region. The back flow has been reduced with forward 
blade. In addition, the effects of tapered hub and diffuser angle 
are combined to minimize the back flow. Figure 14 shows the 
efficiency curves for these design variations. The inducer 
addition clearly shows substantial improvement in the 
hydrodynamic efficiency, at the same time provides a sufficient 
pressure rise to suppress cavitation. 



Figure 14. Hydrodynamic efficiency distribution along impeller 
blade height of various designs. 


Besides improving the pumping efficiency, the design of the 
VAD requires good wall washing near the solid wall by reducing 
the stagnation regions. One of the critical regions for potential 
blood clotting is near the bearing area between rotating and non- 
rotating components. Clotting can be caused in the hub area due 
to either high shear or stagnation depending on the gap and 
configuration of the area. 

Figure 15 shows velocity vectors colored by velocity magnitude 
; for different bearing designs. Design 1 is the original baseline 
design with the cavity width of b. This design showed very high 
shear stresses near the rotating hub face and very stagnant fluid 
region in the lower portion of the cavity. Increasing the cavity 
width gradually up to 8b, circulation in the cavity is increased 
substantially. In order to eliminate stagnant areas in the lower 
portion of the cavity, the hub surface was then tapered which 
reduces the cavity height accelerating the flow near the hub 
region. This resulted in stronger wall washing in the cavity. A 
modified version of this configuration has been adopted in the 
current DeBakey VAD design. 

In Figure 16, areas where the VAD design is improved by using 
CFD analysis is summarized. 



Figure 15. Velocity vectors inside the initial (left) and the final 
bearing geometry (right). 



Figure 16. Contribution of CFD analysis to VAD design. 

The performance of the new design is compared to the baseline 
design in Table 1 where the clinical results were obtained by 
BCM. The hemolysis index reported here shows the amount of 
hemoglobin generated by the pump in grams per 1 00 liters. 
Destruction of the red blood cells results in the release of 
hemoglobin. The new design shows a remarkable improvement 
in performance over the baseline design. There is a 22% increase 
in overall pumping efficiency. The performance of the new 
design was sufficiently improved resulting in a device 
implantable in human. The first human implantation was 
performed successfully in 1998 followed by applications to many 
other patients ever since. 
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Table 1 . Clinical results for baseline and new designs of VAD. 

SUMMARY 

An efficient and robust solution procedure for 3-D turbopump 
analyses and its spin-off application to VAD impeller have been 
presented. The technique solves the viscous incompressible 
Navier-Stokes equations using the pseudo-compressibility 
method. A high-order accurate upwind differencing and an 
efficient implicit scheme are used in the flow solver. The flow 
simulation procedure was developed to analyze an advanced 
turbopump impeller and an SSME configuration. Validated 
solution procedure was then applied to the development of 
DeBakey VAD. Various design improvements were made through 
the use this computational tool. For example, an inducer addition 
dramatically increased pumping efficiency thereby reducing the 
hemolysis to an acceptable level for human use and an optimum 
cavity redesign practically removed thrombus formation in the 
bearing area. Overall the VAD development was expedited with 
the help of CFD technology originally developed for rocket 
pump, thus enabling human implantation. Overall performance 
has been demonstrated through successful human implantations. 
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